/*!
 ***************************************************************************
 * \file   oops.cpp
 * This code solves the heat diffusion equation in 2D.
 ***************************************************************************
 *  \author    Luis M. de la Cruz Salas [ Mon Apr 14 09:28:22 BST 2008 ]
 ***************************************************************************
 * Heat diffusion equation in two dimensions: 
 * \f$ \frac{\partial T}{\partial t} = \nabla^2 T\f$. 
 * \par
 * This equation is solved in a unit square interval \f$ x, y \in [0,1] \f$. 
 * The initial condition is "OOPS"(Object-Oriented Programming for Science) 
 * inside the square. The boundary conditions are \f$ T = 0 \f$ in all walls 
 * of the square.
 *  \par Compiling and running
 *  Modify the variables BASE and BLITZ in the file \c tuna-cfd-rules.in 
 *  according to your installation and then type the next commands:
   \verbatim
   % make
   % ./oops \endverbatim
 * The results of this example are stored in \c DataOOPS.
 * Next figures show these results for the initial and an advanced
 * time of the simulation. The graphics were plotted using OpenDX 
 * with the \c vistempOOPS.net script.
 * \image html  oops1.png "Initial temperature distribution" width=10cm 
 * \image latex  oops1.eps "Initial temperature distribution" width=10cm 
 * \image html  oops2.png "Temperature distribution after some time steps" width=10cm 
 * \image latex  oops2.eps "Temperature distribution after some time steps" width=10cm 
 ***************************************************************************
 */
#include <iostream>
#include <sstream>
#include <cmath>
#include "Meshes/Uniform.hpp"
#include "Equations/ScalarEquation.hpp"
#include "Schemes/CDS.hpp"
#include "Solvers/TDMA.hpp"
#include "Utils/inout.hpp" 
#include "Utils/num_utils.hpp"
#include "Utils/GNUplot.hpp"

using namespace Tuna;

// Scalar field in two dimension
typedef TunaArray<double,2>::huge ScalarField2D;

void initialConditionForTemperature(ScalarField2D &T)
{
    double data[1600];
    ifstream fp("./DataOOPS/OOPS_Data");
    for(int j = 19; j >= 0; --j)
	for(int i = 0; i < 80; ++i)
	    fp >> T(i,j);
}

int main() 
{
    double length_x = 4, length_y = 1;
    int num_nodes_x = 80, num_nodes_y = 20;
    double tolerance = 1e-5, dt = 1.0e-4;
    int frequency = 1;

    StructuredMesh<Uniform<double, 2> > mesh(length_x, num_nodes_x,
					     length_y, num_nodes_y);
    double dx = mesh.getDelta(X);
    double dy = mesh.getDelta(Y);
    mesh.print();

    ScalarField2D T(mesh.getExtentVolumes()); 
    ScalarField2D Tn(mesh.getExtentNodes());
    initialConditionForTemperature(Tn);


    NumUtils::interpolateToVolumes(T, Tn);

#ifdef WITH_DX
    InOut::writeToFile_DX(Tn, 0, "./DataOOPS/temp.", dx, dy);
#endif

#ifdef WITH_GNUPLOT
    InOut::writeToFile(Tn, 0, "./DataOOPS/temp.", dx, dy);
    GNUplot plotter;
    plotter("set contour");
    plotter("set cntrparam level 20");
    plotter("set grid");
    plotter("splot [][][0:1]\"./DataOOPS/temp.00000\" title \"OOPS!\" w l");
    std::cout << "\n >---> Press <1> and then <enter> to continue \n";
    int pause; std::cin >> pause;
#endif
    SparseMatrix< Diagonal <double, 2> > A(num_nodes_x, num_nodes_y); 

    ScalarField2D               b(num_nodes_x, num_nodes_y);

    ScalarEquation< CDS<double, 2> > oops_eq(T, A, b, mesh.getDeltas());
    oops_eq.setDeltaTime(dt);
    oops_eq.setGamma(1.0);
    oops_eq.setDirichlet(LEFT_WALL);
    oops_eq.setDirichlet(RIGHT_WALL);
    oops_eq.setDirichlet(BOTTOM_WALL);
    oops_eq.setDirichlet(TOP_WALL);
    oops_eq.print();

    ofstream error_fp("./DataOOPS/error"), residual_fp("./DataOOPS/residual");
    int iteration = 0;
    double errorcillo = 1.0, residual = 1.0;
    
    while (errorcillo > tolerance) {
	oops_eq.calcCoefficients();
       	Solver::TDMA2DX(oops_eq, tolerance, 10);
	errorcillo = NumUtils::calcErrorL1(T, oops_eq.phi);
	residual = oops_eq.calcResidual();

	error_fp << iteration << "\t" << errorcillo << endl;
	residual_fp << iteration << "\t" << residual << endl;
	
     	cout << "\n ---- Iteration " << ++iteration
	     << " Residual = " << residual 
	     << " Error = " << errorcillo;
	
	oops_eq.update();
	if ( !(iteration % frequency) ) {
	  NumUtils::interpolateToNodes(Tn, T);

#ifdef WITH_DX
	  InOut::writeToFile_DX(Tn, iteration, "./DataOOPS/temp.", dx, dy);
#endif

#ifdef WITH_GNUPLOT
	  InOut::writeToFile(Tn, iteration, "./DataOOPS/temp.", dx, dy);
#endif
	}
    }
    
    
#ifdef WITH_GNUPLOT
    std::cout << "\n >---> Press <1> and then <enter> to plot the solution \n";
    std::cin >> pause;	
    for(int i = 0; i < iteration; i += frequency) {
      std::string name = "./DataOOPS/temp.";
      std::ostringstream extnumber;
      extnumber.width(5);
      extnumber.fill('0');
      extnumber << i;
      name += extnumber.str();
      std::string command = "splot [][][0:1] \"" + name + "\" w l";
      plotter(command);
    }
#endif

#ifdef WITH_DX
  std::cout << "\n\n NOTE: The data were written in the [./DataOOPS] directory in DX format";
#endif

    return 0;
}




